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We consider the question of the amorphization of metallic alloys by melt quenching, as predicted by molecular 
dynamics simulations with semi-empirical potentials. The parametrization of the potentials is discussed on 
the example of the ternary Cu-Ti-Zr transition metals alloy, using as reference the ab-initio simulation. 
The pair structure in the amorphous state is computed from a potential of the Stillinger Weber form. The 
transferability of the parameters during the quench is investigated using two parametrizations: from solid 
state data, as usual, and from a new parametrization on the liquid structure. When the adjustment is made 
on the pair structure of the liquid, a satisfactory transferability is found between the pure components and 
their alloys. The liquid structure predicted in this way agrees well with experiment, in contrast with the one 
obtained using the adjustment on the solid. The final structure, after quenches down to the amorphous state, 
determined with the new set of parameters is shown to be very close to the ab-initio one, the latter being 
in excellent agreement with recent X-rays diffraction experiments. The corresponding critical temperature of 
the glass transition is estimated from the behavior of the heat capacity. Discussion of the consistency between 
the structures predicted using semi-empirical potentials and ab-initio simulation, and comparison of different 
experimental data underlines the question of the dependence of the final structure on the thermodynamic 
path followed to reach the amorphous state. 


I. INTRODUCTION 

Since their discovery in the late eighties^, multicom¬ 
ponent bulk metallic glasses (BMG) have been the sub¬ 
ject of numerous experimental and theoretical studies. 
The unique physical, mechanical, and corrosion proper¬ 
ties of these novel materials indeed enable a variety of 
applications^Ti^. A particular attention has been paid to 
the BMGs based on late transition metals, especially in 
the Gu-based family^i^, including alloys with Zr—Ti^, 
or AliS. These studies have shown that the glass forming 
ability (GFA) of an alloy and its properties in the amor¬ 
phous state depend on several factors such as the nature 
of its components, its composition, or the cooling rate. 
Simulations are then helpful for a systematic exploration 
of the parameters space, including domains that are dif¬ 
ficult to study experimentally. They are also useful for 
understanding the underlying physics at an atomic scale. 
Most often, they consist in classical molecular dynam¬ 
ics (MD) with semi-empirical atomic potentials, such as 
the embedded atom model (EAM)ii, the tight binding, 
second moment approximation (TB-SMA)^^, and Finnis- 
Sinclair (F-S) potentials'^ and their variants. Their pa¬ 
rameters are adjusted from available structural and ther¬ 
modynamic data. EAM potentials for fee elements are 
given for example in Ref. Q and n-body potentials for 
ternary alloys are reviewed in Ref. [la. For the Gu family 
see Ref. [l6| for Gu-Zr-Ag, or Ref. Il7l for Gu-Zr alloys. 
For the latter, the parameters were determined from a 
combination of ab-initio calculations and experimental 
data. Other methods, such as the force matching tech¬ 
nique in which ab-initio potential energies are used as 
input, have also been used (see for example Refs. [T^ for 


the Gu-Zr-A1 alloy and fl^ for Gu-Zr). Garefully adjusted 
potentials can then be used for practical applications^^, 
like searching optimized compositions^^. 

The central question in this classical route remains the 
transferability of the parametrized force fields, to differ¬ 
ent state points or compositions, for example. The trans¬ 
ferability, especially for the n-body potentials has been 
discussed e.g. in Ref. [l^for fee metals in the liquid phase 
(TB-SMA), Ref. lU for Ni-Al (EAM), or Ref. ^ for 
liquid Pd-Ni alloys (modified Sutton-Ghen^^ potential). 
Furthermore, the data used for the fits - say the cohe¬ 
sive energy, lattice constants, etc. - are usually relative 
to the solid phase. It is then not obvious that the same 
parameters will also describe the properties of liquid and 
amorphous metals^, which lack the periodicity which is 
important for the electronic structure of the solid^i. In 
related areas, this question is discussed in Refs. [2^ and 
1^ or [23 for the use of this approach in thermodynamic 
integration methods. 

As an alternative, one may consider first-principles 
simulations for which the question of transferability does 
not arise. Their computational cost however makes them 
often unfeasible without resorting to supercomputers. 
This holds even with non all-electron ones, which use 
pseudopotentials and approximate functionals. This in¬ 
cludes paths involving a wide range of variation of the 
parameters, such as in cooling rate studies^^ - see also 
Refs. [ 2 ^ and [s^- They are thus used occasionally, say 
for supplementing the information required to fit the 
phenomenological potentials2£i^y2i. To this end, when 
the amorphous state is obtained by quenching a liquid, 
namely through a path between two disordered states, it 
may be preferable to adjust the parameters on the liquid 
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properties. 

Since there is no systematic means to devise the op¬ 
timum transferability (for alternatives see for example 
Ref. [2^ . new elements collected in representative cases 
are useful to ascertain this question. In this work, 
we shall illustrate this on the example of the ternary 
Cu6oTi2oZr2o alloy which has been studied recently by 
experiment (see references in Refs. [HL IsHj - ldOh and simula¬ 
tion with parametrized potentials'^— for its importance 
as a model system of transition-metal-based BMGs, and 
its technological relevance. To our knowledge, this al¬ 
loy has not been studied from first-principles simulation, 
in contrast with Cu-based binary alloys (see for example 
Refs, li^andli^for Cu-Zr). Using molecular dynamics at 
the Born-Oppenheimer level (BOMB) and density func¬ 
tional theory (DFT)^, as implemented in the Quantum 
Espresso {QE) package^, we present here such a study 
intended to highlight some questions raised by the use 
of locally adjusted atomic potentials in multicomponent 
glass forming alloys (a recent review of ab-initio molecu¬ 
lar dynamics (AIMD) methods applied to glass formation 
can be found in Ref. li^ . This will be done by com¬ 
paring the pair structure determined from the ab-initio 
and classical routes, in the liquid and at ambient tem¬ 
perature which is well below the critical temperature of 
the glass transition of this alloy. By adjustment to the 
high-temperature structure of the alloy determined from 
AIMD, a new set of parameters of the effective potential 
is determined. Their transferability is investigated by 
comparison with the ab-initio structure at the final tem¬ 
perature. The corresponding critical temperature of the 
glass transition is estimated from the behavior of the heat 
capacity. The stability of the structure predicted using 
the semi-empirical potential with respect to the AIMD 
simulation is discussed. Comparison with other simula¬ 
tions and available experiments is made to underline the 
importance of thermodynamic path followed to reach the 
amorphous state. 

This paper is thus organized as follows: in section II, 
we detail the methodology we use for this purpose. In 
section III, we present the main results concerning the 
ab-initio (g^J’(r)) and classical {gij{r)) radial distribu¬ 
tion functions (rd/s) and the estimated glass transition 
temperature. The system size dependence is finally dis¬ 
cussed. We conclude this paper by a summary of the 
main conclusions. 


II. METHODOLOGY 
A. General method 

The question of transferability is actually independent 
of the reality of the reference data used to parametrize 
the potentials. As we could not find experimental data 
for the pair structure of Cu6oTi2oZr2o at high tempera¬ 
ture, we decided to use as “experimental” data the rd/s 
gfj{r) of the alloy obtained from AIMD in the liquid 


state. We postpone to the next section the question 
of its relationship with the actual structure of the liq¬ 
uid alloy. The “classical” rd/s gij{r) are determined by 
MD simulation using the simplified (i.e. without 3-body 
terms) Stillinger- Weber— (SW) potential used by Han 
and Teichler— : 

(j^swir) = a( - 1 ] exp (- ^ -] , 0 < r < —. 

\{ar)^ J \ar — ai J a 

( 1 ) 

These simulations will be referred to as /)sw-MD. The set 
of parameters of the SW potential adjusted using data 
from the solid^ will be designated as {a®} and those 
from the liquid, determined as {a(}. The parameters 
{a(} are varied until the pair structure determined with 
the ^sw-MD simulation {gij{r)) agrees reasonably well 
with the ab-initio one {gfj(r)). Once this is achieved, 
the alloy is cooled down to T = 300 K. Starting from 
the last equilibrated MD configuration so obtained, the 
AIMD is run at ambient temperature, until satisfactory 
equilibration is reached. This is followed by a series of 
accumulation steps for obtaining the statistical averages. 
A first comparison can thus be made between the two 
final sets of rd/s at T = 300 K, gfj{r) and gij{r). As a 
test of the role of the initial configuration, we perform an 
“instantaneous” quench by using as input at T = 300 K 
a configuration equilibrated in (or near) the liquid state. 
One may then compare the rdf determined with finite 
and infinite cooling rates (with necessary caution in the 
latter case). We also used the last configuration obtained 
with the {a®} parametrization. Before presenting these 
comparisons, we give below some technical details of the 
classical and ab-initio runs. 


B. Computational details 
1. Classical MD 

The extensive ^sw-MD runs were performed using 
the LAMMPS package^ in the isobaric NPT ensem¬ 
ble using the Nose-Hoover integration, with a time step 
dt = 0.0025 ps. Over a total length of about 6 10® steps, 
the last 2 10® were used to compute the averages. All 
potentials were cut at r = a/a and finite size effects were 
investigated by changing the particle number from 260 to 
1372. The quenches were performed in the NPT ensemble 
(at p = 0), either with a stepwise change of temperature 
or a continuous one. Different cooling rates were consid¬ 
ered (besides an instantaneous quench): 3 10^® Ks“^ and 
a 10 times slower one. 


2. Ab-initio MD 

The AIMD simulations were run using of the plane- 
wave self consistent field (PWscf) code in the QE pack¬ 
age with no modification other than adding the compu- 




3 


tation of the six rdfs gfj{r) at regular time intervals (say 
every twelve hours) to monitor equilibration and produc¬ 
tion steps directly on the pair structure. The ions dy¬ 
namics used a time step of dt = 30 a.u. (1.45 10“^ ps), 
the equations of motion being integrated with the Ver- 
let algorithm. Between N = 128 (pure components) and 
N = 240 — 260 (ternary alloy) particle numbers were 
considered, and a simple velocity rescaling method was 
used to fix the temperature. (N,V,T) simulations were 
performed, but, to bracket the desired zero pressure, the 
volume was initially varied starting from the value for 
an ideal mixture (see below). A series of short runs 
(about 200-300 steps) with different box sizes were made 
to monitor the fluctuation of the pressure above and be¬ 
low p = 0. Full NPT simulations using the variable 
cell-shape method implemented in the QE distribution 
indeed proved too costly on the rather small 24-core ma¬ 
chine we used. 

The accuracy of the AIMD simulation depends mostly 
on the quality of the electronic energy computed in the scf 
cycles. The quality of the exchange-correlation functional 
used in the DFT treatment of the electronic structure and 
some parameters such as the number of plane waves used 
in the expansion of the Kohn-Sham orbitals are known to 
be important, but they were found not to be critical here 
(i.e. for the determination of the pair structure). After 
several trials, an energy cutoff of 30 Ry and a density 
cutoff of 240 Ry, in the range of the recommended val¬ 
ues, were found sufficient to ensure a satisfactory conver¬ 
gence. Ultra-soft (USP)^ pseudopotentials from the QE 
library^, all with the PBE functional of Perdew et alM- 
were used (with Z = 11,12,12 valence electrons for Cu,Ti 
and Zr). While more accurate^ projector augmented- 
wave (PAW)^ pseudopotentials are available, our tests 
did not indicate a significant effect at the level of the rdfs 
(see below the discussion for pure Cu). The convergence 
threshold for the scf cycles was 10“® Ry for typical en¬ 
ergies of 3 10^'’ Ry. Due to the large supercell size used 
(typically 15 A), calculations were performed at the P 
point only (except for the density of states discussed in 
section ITlI Cl) . In these conditions and using standard val¬ 
ues (for metals) for the other settings of the PW code, 
one scf cycle for a total number of 2964 electrons (for 260 
atoms) takes on our machine roughly 2300 s. 

III. RESULTS 

A. High temperature structure 
1. Liquid copper 

To test various parameters in the AIMD runs, we com¬ 
puted the rdf g’cui'’’) Pill's Cu. N = 128 parti¬ 
cles were placed in a cubic box with a supercell size 
L = 23.3 a.u. to achieve the experimental equilibrium 
density p = 0.074 at atmospheric pressure^, g'cuir) 
was computed in the liquid at T = 1623 K to allow com¬ 


parison with the AIMD of Ganesh and Widom^ (PAW 
pseudopotentials) and the experimental data of Ref. 

The AIMD parameters are those indicated above. The 
result is shown in figure O- 



FIG. 1. Radial distribution function of liquid Cu at T = 
1623 K. 

Symbols: experimental X-ray diffraction data^ interpolated 
to 1623 K. Full line: simulation of Ref. (PAW). Dots: this 
work (USP). 

To test the sensitivity to the pseudopotential, both 
USP and PAW ones were used. The value of the pressure 
is actually sensitive to the pseudopotential, but the dif¬ 
ference between USP and PAW pseudopotentials (a few 
kilobars here) has a negligible effect on the equilibrium 
density and hence on the rdf. With the USP, the peak is 
higher by two percents but the oscillations are of similar 
quality as those with the PAW potentials. The simi¬ 
lar agreement between USP and PAW calculations with 
experiment led us to use the USP since the scf cycles 
are then less time-consuming. Computing an accurate 
rdf for Cu has actually no novelty per se since this can 
be achieved using a much faster parametrized version of 
the generalized pseudopotential theory for the noble met¬ 
als, as known since the mid-eighties^. This calculation 
was only intended to check various parameters in the 
scf cycles, besides giving an indication of the expected 
structure in the alloy. A previous study by Jakse and 
Pasture R^i^^ indeed showed a relative insensitivity of the 
ab-initio Cu-Cu rdf (at least in the first peak) to copper 
concentration in the binary Cu^Zri-j, alloys. We could 
not find a similar study for Cu-Ti but the data in Refs. 
1^ and do not indicate qualitative effects of alloying 
with Ti. The important point then is that since Cu is 
the majority component in Cu 6 oTi 2 oZr 2 o, the rdf of Cu 
should be similar in the alloy and in the pure compo¬ 
nent, which is very well described by AIMD. This is also 
a good indication that the structure determined from the 
simulation of a small periodic system is representative of 






4 


the one in the true alloy, as discussed later. 

2 . Liquid alloy 

The partial rd/s of Cu 6 oTi 2 oZr 2 o at 1600 K deter¬ 
mined from the AIMD with ultra-soft pseudopotentials 
are shown in figure They were obtained with N — 
260 particles (156 -1- 52 -1- 52 atoms for Cu, Ti, and Zr 
respectively) to keep the CPU time in a tolerable range. 
The system was initially prepared as follows: an initial 
density was computed from the known experimental den¬ 
sities of the pure liquid components, assuming an ideal 
mixture behavior. This gives p ^ 0.062 A“^. For the 
corresponding volume, a cubic box is filled with the N 
particles randomly distributed on the vertices of a cfc 
lattice (plus a small random displacement), to form the 
supercell. Its size is then varied about the initial value 
L ~ 30.5 a.u. to bracket the zero pressure at 1600 K. 
Within a tolerance of a few kilobars, L converged rather 
rapidly, the final value being L = 30.914 a.u. The cor¬ 
responding volume was then used to continue the AIMD 
runs in the NVT condition. After about 500 steps to 
monitor equilibration, 1500 steps of accumulation gener¬ 
ate, in about 45 days of computation, reasonably smooth 
curves for gcu-Cu, as well as for the cross ones involv¬ 
ing Cu, somewhat less so for the minority species Ti 
and Zr . Note for the latter a significantly larger effec¬ 
tive diameter and a “smoother” g{r), clearly calling for 
a refined parametrization, possibly with the three-body 
terms (with the present model, we recently reached a to¬ 
tal of 2475 steps which showed no significant effect aside 
from a progressive smoothing of the data). As a check of 
equilibrium, we noted that the total energy is very sta¬ 
ble, with a very small relative fluctuation and the pres¬ 
sure recorded during the last 839 steps fluctuates about 
a slightly negative but well-defined average valued. 

The set of parameters {a(} was then determined by 
(hand) adjustment on these ab-initio rd/s, in the liquid 
ternary alloy. They are given in table I. 



A 

(eV) 

a 

^(A) 

ai 

n 

Cu-Cu 

1.135 

(0.485) 

2.290 

(2.275) 

1.681 

(1.681) 

9(9) 

Cu-Ti 

1.725 

(1.695) 

2.325 

(2.300) 

1.805 

(1.794) 

7(7) 

Ti-Ti 

1.364 

(1.588) 

2.380 

(2.350) 

1.960 

(2.056) 

9(4) 

Cu-Zr 

1.942 

(1.943) 

2.450 

(2.496) 

1.792 

(1.792) 

8(8) 

Ti-Zr 

1.959 

(2.722) 

2.530 

(2.481) 

1.900 

(1.968) 

7(3) 

Zr-Zr 

1.695 

(3.655) 

2.710 

(2.646) 

1.950 

(1.855) 

9(3) 


TABLE I. Parameters {a(} of the SW potential adjusted on 
the ab-initio liquid structure. The values between parentheses 
shown for comparison are the parameters {al} adjusted on the 
solid, from Table II of Ref. [dj. 


The corresponding g^j are drawn as lines in figure (P|) . 



FIG. 2. Radial distributions functions in the Cu 6 oTi 2 oZr 2 o 
alloy at T = 1600 K. 

Symbols: AIMD; lines: /sw-MD (with parameters {a\} as in 
table I). 


It is stressed that this set of parameters is by no 
means unique. The goal was not to determine the best 
parametrization but to investigate whether one that gives 
a reasonably accurate initial structure (in comparison to 
a given reference) keeps its quality after a quench down 
to the amorphous state. For this reason, the much longer 
runs that would have been needed to improve the statis¬ 
tics for the minority species were not attempted. With 
this caveat in mind, it is clear from the table that the val¬ 
ues determined from the adjustment to the liq uid struc¬ 
ture differ considerably from those of Ref. We note 
in particular the large difference for the exponent pa¬ 
rameters which affect the short-range behavior of such 
effective potentials (see above the remark for Zr). This is 
a clear evidence of non-transferability from the solid to 
the liquid phase. 
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3. Transferability of the parametrization {a\} at high 
temperature 

The transferability of the parameters {a*} in the liquid 
state was first checked on the pure components. We show 
in figure ([3]) the rd/s for the pure components computed 
with the parameters {a*} adjusted on the ab-initio rd/s 
in the alloy, along with experiment for pure Cu^, Zr— , 
and Ti^. 



r [nm] 


FIG. 3. Radial distribution functions of liquid Cu (1623 K), 
Ti (1973 K) and Zr (2173 K). 

Symbols: experiment (Cu and Zr: Ref. Is^ . Ti: Ref. Ig^I : full 
curves: /sw-MD, {a*}; dotted curve: /sw-MD, {a|} for Cu; 
dashed lines: AIMD (from Ref. for Zr at 2500 K). 

Similarly to the ab-initio ones, the rd/s for the pure 
components determined with the parameters {a*} ad¬ 
justed on the alloy are in very good agreement with ex¬ 
periment. On the contrary, with the parameters {a®} de¬ 
termined from lattice constants and cohesive energy data 
as in Ref. Idll g{r) differs significantly from experiment 
(see figure ([3]) for Cu). The situation is thus different 
from that of the pure components since it was observed 
in Ref. [13 for example that tight-binding potentials give 
a reasonable description of the dynamic properties of sev¬ 
eral liquid metals, in spite of having been parametrized 
on the basis of solid-state data. Determining fully trans¬ 
ferable parameters from the solid to the liquid and vice- 
versa remains thus an important task in the future. 

From the good transferability between the ternary al¬ 
loy and the pure components, one expects a similar qual¬ 
ity of {aj} for the binary alloys. This was checked for 
the Cu-Zr alloy for which data for the pair structure are 
available. In figure 0 we show the comparison with 
the ab-initio results of Jakse and Pasturel^. It clearly 
confirms the expectation. We did not investigate other 
binary alloys formed with Cu, Ti and Zr, but we see a 
priori no reason why the parametrization should then be¬ 
have differently. 



r [nm] 

FIG. 4. Radial distribution functions in liquid Cu 8 oZr 2 o at 
T = 1500 K. 

Symbols: AIMD of Ref. l45l . Full lines: /sw-MD, {a\}. 

A test of different nature is to compare the predicted 
structure dir ectly with experiment. We did this with the 
data of Ref. l65l for the Zr 7 oCu 3 o alloy at T = 1453 K. 
This constitutes a severe test since Cu is a minority com¬ 
ponent for this composition, gtot then has a large con¬ 
tribution from gzr-Zr and gzr-Cu- The result is shown 
in figure 0. The agreement between experiment and 
simulation is quite satisfactory, especially in view of the 
fact that the adjustment of the parameters for Zr was 
made with only 52 atoms in the ternary alloy. Our re¬ 
sults are comparable to those of the Monte Carlo simula¬ 
tion of Harvey et al^ using a modified EAM formalism. 
Studying, possibly with refined parameters, the amor¬ 
phous structure of binary Cu-Zr alloys (like the shoulders 
in the first peak of gtotr^, its sensitivity to composition^, 
or the effect of the cooling rate^) is left for future work. 

To conclude this section, these comparisons with ex¬ 
periments and other simulations show that when the ad¬ 
justment is made on the structure in the liquid state, the 
transferability from the alloy to the pure components and 
the binary alloys is very good. Thus we have the desired 
starting point for investigating whether it remains so, af¬ 
ter a quench down to ambient temperature. 


B. Amorphous alloy 

1. Structure from classical and AIMD simulations 

Using the parameters {a\} in table I in the ^sw-MD 
simulation, we perform a quench from T = 1700 K to 
T = 300 K, first at the cooling rate of 3 10^° Ks“^ used 
in Ref. Ull The rd/s at the final temperature are shown 
in figure ®. We first observe the split second peak. 
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FIG. 5. Total radial distribution functions in the liquid 
ZrroCuso alloy at T = 1453 K. 

Solid circles: experimenls^; full curve: 0sw-MD, parameters 
{a*}; dotted line: same without the volume term in the pres¬ 
sure to test sensitivity to the average density. 


typical of the amorphous solid. The first peak is also 
much higher than in the liquid and also more pronounced 
than the one predicted with the parameters {af} (figure 
(6) in Ref. SI. 

To test the resulting structure at T = 300 K, we de¬ 
termined the AIMD gfj{r) starting from the final con¬ 
figuration of the (()sw“MD run (below designated as con¬ 
figuration (1)). To reduce CPU time, as indicated, NVT 
simulation was used, for the average volume of the (j)sw- 
MD run at zero pressure. The corresponding cell size 
is 29.837 a.u. After runs of length similar to the high- 
temperature ones (about 500 -I- 1500 steps), sufficiently 
smooth curves were obtained. We first observe in fig¬ 
ure ([6]) that the rd/s involving Cu are very close in the 
ab-initio and classical MD simulations. The larger dis¬ 
crepancy for Ti-Ti and Zr-Zr is mostly a consequence of 
the smaller statistics in the ab-initio data, obtained with 
only 52 atoms. This is confirmed by the smooth behav¬ 
ior of the partial rd/s involving Cu in figure, which are 
averaged over a larger number of pairs. 

To check equilibration, the pressure was recorded dur¬ 
ing the last 750 steps. It oscillates about a stable (albeit 
slightly too high) average value2£. This should not af¬ 
fect much the pair structure, since the amorphous al¬ 
loy is rather dense at T = 300 K {p = 0.0661 
to be compared to p = 0.0594 at T = 1600 K). 
Since the dynamics is particularly slow in the amorphous 
state, sensitivity to the initial configuration was tested 
by starting the AIMD run from two other equilibrium 
configurations : (2) the AIMD one at T = 1600 K, after 
2000 steps - this amounts to performing an instantaneous 
quench (the positions were rescaled by the density ratio 
(0.059/0.066)1/3); (3) the (/sw-MD one at T = 300 K, 



FIG. 6. Partial rdfs in amorphous Cu 6 oTi 2 oZr 2 o at 300 K. 
Symbols: AIMD; full lines : /sw-MD, {a*}; dashed line 
for Cu: MD simulation^ with tight-binding potentials for 
Cu5oTi25Zr25. 
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after 610® steps with the parameters {af}. The latter 
differs from configuration (1) (the one with the parame¬ 
ters {a*})- The results shown in figure ([7]) indicate that 
(i) the initial condition plays a role, which is understand¬ 
able for amorphous states, but the final AIMD rd/s are 
not strongly affected, (ii) the ab-initio simulation is not 
frozen in the initial condition as evidenced by the com¬ 
parison between the final AIMD rd/s at 300 K and those 
corresponding to the initial conditions: g^„(r) is very dif¬ 
ferent from that in the liquid (configuration (2)), even for 
a liquid having the density of the amorphous alloy2i; it 
is also clearly different from that in the amorphous state 
as predicted by (/)sw-MD with parameters {of} (configu¬ 
ration (2)) - see also inset in figure ([5]) for gtot- Another 
indication is the mean squared displacement. The AIMD 
value estimated on the (limited) number of available con¬ 
figurations, Ar^ ~ (0.06 ± 0.01) A^, is compatible with 
the one determined by classical MD^^., both being differ¬ 
ent from the one obtained by classical MD with parame¬ 
ters {af}. This suggests that the latter parametrization 
leads to a faster dynamics in the frozen state. A larger 
diffusivity is also clearly evidenced by the behaviour of 
Ar^(f) in the liquid, as determined by classical simula¬ 
tion over a wider time range (the same should hold also 
in the undercooled states)^. 

The behaviors of the AIMD simulation with configu¬ 
rations (2) and (3) are then clearly distinct from the one 
observed with configuration (1), that leads to the agree¬ 
ment with experiment. While one might naturally fear 
freezing for relatively short AIMD runs, this dependence 
on the initial configuration shows that the microstates 
explored by the AIMD simulation remain close to those 
of the classical MD only when they correspond to the cor¬ 
rect free energy minimum^i. This is reassuring in view 
of the agreement between the classical and ab-initio rd/s 
in figure (|5]) when using the parameters {a(}. 

Regardless of the connection between the structure 
predicted by simulation and the one in the macroscopic 
alloy (see figure ®), one first important result is thus 
the good transferability of the potentials from the liquid 
to the amorphous state. 

A second important observation is the noteworthy sim¬ 
ilarity of the structure predicted for Cu 6 oTi 2 oZr 2 o to that 
of the Cu 5 oTi 25 Zr 25 alloy, with similar composition, stud¬ 
ied by Senturk Dalgic and Celtek^ using TB potentials 
in the simulation (dashed line in figure (IHl)). The amor¬ 
phous state of this alloy was also obtained from a quench 
of the high-temperature liquid (the predicted densities at 
300 K, p ^ 66.13 nm“® and p ~ 63.09 nm“®, respectively, 
are close as well). Such consistent predictions for simi¬ 
lar systems deduced from different potentials, adjusted 
on different data, and the fair agreement with the AIMD 
for Cu 6 oTi 2 oZr 2 o constitute a strong evidence that the 
predicted structure is the actual one, when the amor¬ 
phous state is reached through a rapid quench from the 
liquid. 

This is confirmed by comparison with experiment in 
figure ([5]). The agreement between simulation and the 



FIG. 7. Total ab-initio radial distribution at T = 300 K from 
different initial configurations. 

Full curve: last configuration with {a*}; open diamonds: 
last configuration with {af}; solid diamonds: “instantaneous” 
quench. 

recent X-rays diffraction data of Durisin et aZ4S is nearly 
perfect. This is quite remarkable, considering the absence 
of adjustable parameters in ab-initio simulations (aside 
from the use of pseudopotentials, which are purely atomic 
properties). 



r [nm] 

FIG. 8. Experimental and simulated total rdf of amorphous 
Gu 6 oTi 2 oZr 2 o at T = 300 K. 

Solid circles: XRD data of Durisin et al.—, crosses: AIMD; 
full curve: /sw-MD. 

This pair structure is clearly different from the one pre¬ 
dicted by Han and Teichler using the {of} parametriza¬ 
tion of the SW potential (figure ([5])). Since these pa¬ 
rameters fitted on the ordered solid cannot describe the 
















FIG. 9. Effect of the parametrization on the total rdf at 
T = 300 K. 

Dotted line: i^sw-MD with {af}, full curve: same with |al|. 
Squares: unpublished XRD data of Ref. shown in Ref. Idll . 


structure of the liquid either, as shown above, both the 
high- and low-temperature limits suggest a possible prob¬ 
lem with this parametrization. The puzzling point how¬ 
ever is that the same parametrization predicts a total rdf 
gtot{T) in rather good agreement with the (unpublished) 
experimental data^^ shown in Ref. (figure ([^ ). 

Given the excellent agreement between simulation and 
the experimental data of Durisin et a/JS, which are 
clearly different from those of Ref. and excluding 
artifacts, one possibility is that the experimental rd/s 
in figure (|5]) are different because they actually corre¬ 
spond to different amorphous states. This idea of some 
dependence on the thermodynamic path through which 
the amorphous state is reached is found for example in 
Refs. 1^ [76l - [^ The parametrization with {af} - which 
is not validated by AIMD - might thus correspond to 
a metastable state close to the one reached in Ref. 

(the path for these data^^ is likely the same as in Ref. 

which does not show the rdf). A slight change in 
the thermodynamic parameters or in the parametriza¬ 
tion should thus drive the system away from a local free 
energy minimum. To test this, we first performed a 10 
times slower quench (3 10® Ks“^), with the parametriza¬ 
tion {a(}. We indeed find a detectable effect on gtotir). 
This differs from the study in Ref. of the amorphous 
Cu 46 Zr 54 metallic glass, which showed no significant ef¬ 
fect of the quenching rate in the higher range investi¬ 
gated. The authors used TB-SMA n-body potentials 
parametrized according to the method of Ref. in 
which the force-field parameters were obtained from a 
fit to first-principles/DFT calculation, as we do it here. 
The slight effect of the cooling rate we found is however 
insufficient to explain the difference between the two ex¬ 
periments. A path dependence is also suggested by the 


solid solution model in Ref. It is difficult to compare 
our results with those for the Cu 7 oZra,Tii_a; systems, but 
the great sensitivity to the concentration of the minority 
species shown clearly testifies for the importance of the 
path leading to the amorphization - the solid solution re¬ 
maining crystalline at some concentrations, for example 
(see the discussion of figure (8) in Ref. Id^ . 

In the numerical quench experiment, another source for 
the observed discrepancy is the initial structure, which, 
as discussed above, is quite different with {of} and {a(}. 
As a test of the sensitivity to the parametrization, we 
relaxed the adjustment of the strength parameter for Cu 
to start with a structure that is in-between the “correct” 
one (i.e. with {a*}) and the one obtained with {a®} (re¬ 
call that the latter predicts a less structured liquid). As 
expected, the resulting gtot (?') at 300 K (not shown here) 
is indeed in between the one obtained with the initial 
{a(} parameters (or experiment) and the one with {a®}. 
Another indication suggesting metastability of the state 
described by the data of Ref. is the fact that the 
AIMD run quickly departs from configuration (3) (see 
also figure ©)• More experiments would thus be useful 
to ascertain the structure in amorphous states reached 
through different thermodynamic paths. 

2. Glass transition temperature 

The simulation of a small periodic systems is not ex¬ 
pected to give an accurate description of the amorphiza¬ 
tion of the real ternary alloy, besides other questions such 
as surface effects in the experiment. Its computational 
cost for the required larger systems discouraged us from 
attempting it in this work. Nevertheless, to test the idea 
of adjusting the parameters in the liquid, we used the 
classical (/sw-MD route with the parameters {a\} to es¬ 
timate the critical temperature of the glass transition. 
This was done from the behavior of the constant-pressure 
heat capacity Cp, shown in figure (1101) . 

We thus find Tg ~ 787 K, slightly above the value de¬ 
termined with the parameters {a®} {Tg ~ 750 K). This 
is consistent with the diffusivity being lower with the pa¬ 
rameters {a(}. Both are above experimental estimates 
~ 710 K in Ref. i or ~ 695 K in Ref. 

. Comparison with experiment is however not imme¬ 
diate for several reasons, including the cooling rates that 
are achievable in simulation^. The higher cooling rates 
in the simulation should indeed overestimate the critical 
temperature. A clarification of this aspect on the basis 
of ab-initio simulation is left for future work. 


C. Finite-size effects 

The main question we investigated here is the trans¬ 
ferability at different temperatures and compositions of 
parametrized potentials for a representative ternary alloy, 
using, for computational convenience, a rather small sys- 
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FIG. 10. Heat capacity versus temperature. 

Full curve: (jisw-MD with {a*}; dashed curve: same with 

{af}. 

tern. This question is actually independent of the close¬ 
ness of the properties of the small periodic system to 
those of the macroscopic alloy. Aside from the question 
of statistics, using a small system does not indeed con¬ 
stitute a limitation as long as the goal remains to test 
transferability. To have an idea of the influence of the 
system size, some calculations were nevertheless repeated 
using N = 240 atoms, and no significant change of the 
pair structure was found. 

Actually, the possible impact of using a small periodic 
system has two different aspects: the first one concerns 
the variation of the structure with the number of atoms 
considered. The second one is the metallic character of 
a small system, especially when the minority species are 
represented by a few tens of atoms. 

The first aspect can be important due to the lack of 
long-range order of amorphous materials which requires 
using blocks containing at least several hundreds of par¬ 
ticles (see Ref. for Cu-Ti). To investigate this, we 
repeated the classical simulations with much larger sys¬ 
tems. As we found no significant effect on the rdf, the 
artificial periodicity related to a small system size should 
not be a qualitative limitation, as far as the average struc¬ 
ture is concerned. This is confirmed by the very good 
agreement with the recent X-ray data shown in figure 
([5]). On the other hand, the ab-initio rd/s were obtained 
in NVT conditions for a rather small system having the 
average density at zero pressure in the classical MD. This 
gives a slightly too high ab-initio pressure while the small 
system size makes it slightly anisotropic. But due to the 
very steep variation of the pressure with density only 
very accurate values for the latter might require larger 
systems. 

Other structural properties might however be much 
more sensitive to the size of the system (see e.g. Refs. 
HI, [sO, and HO). In Ref. [s^, for example, it has been 


noticed that selective minor additions can dramatically 
improve the GFA of binary metallic glasses. In particu¬ 
lar, the effect of composition on the glass formation of the 
Cu-Ti-Zr alloys has been discussed in Refs. Isbl - I^ and 
1^ For the family of Cu-based alloys, a recent study^ 
based on the parametrization of Ref. [13 evidenced the 
subtle structural evolution that occurs during cooling. 
They pointed out that the presence of icosahedrally co¬ 
ordinated clusters and their tendency to form networks is 
insufficient to explain glass formation at all compositions 
in the Cu-Zr binary system. Besides the Cu-Ti-Zr al¬ 
loy considered here, there are several studies of Cu-based 
ternary alloys, such as Cu-Zr- AliS and Cu-Zr- Agi^. In 
Ref. Ra a many-body potential was developed using the 
embedded atom method (EAM) on the basis of ab-initio 
calculations^^. They pointed out the coupling between 
chemical and dynamical heterogeneities, which appears 
to play a crucial role in the improved GFA of this alloy 
and the Cu family of alloys studied in Ref. [s^- Studying 
such size-dependent structural effects by ab-initio sim¬ 
ulations would likely require larger samples than those 
considered here. 

The second question relates to how far the electronic 
structure would have been different had we used a larger 
system, due to the actual sensitivity of the metallic char¬ 
acter to the number of atoms in a small sample, as dis¬ 
cussed in the literature on the related field of metallic 
clusters (see e.g. Refs. and references therein). 

For the actual size used here, an idea can be formed from 
the system size dependence of typical electronic proper¬ 
ties such as the density of states (DOS), the Fermi level, 
or the energy per particle. Figure m shows the instan¬ 
taneous DOS in equilibrated configurations of the nuclei 
at 300 K: two for N = 240 atoms and one for N = 260. 
In the three cases, the DOSs are very similar, being dom¬ 
inated by the d-band-like contribution of Cu in the range 
—6eV ^ e — ep ^ —2 eV. The corresponding Fermi en¬ 
ergies are ep = 13.4181, 13.4228, and 13.3775 eV, and 
the total energies per particle are E/N = —1578.922, 
— 1578.919, and —1578.887 eV, respectively. A quite sim¬ 
ilar behavior for Cu-Zr binaries is shown in figure (4) of 
Ref. HU This shows that the variation with system size 
is of the same amplitude as the variation between two 
configurations having the same number of particles. Av¬ 
eraging over a large number of configurations will smooth 
even more the tiny differences found here (see also Ref. 
Is^ . Thus the systems we used seem large enough also 
for the electronic structure at fixed configuration of the 
nuclei. 


IV. CONCLUSION 

We have discussed in this work the transferability of a 
simplified model potential of the Stillinger-Weber form, 
between the pure components and the alloy in the liq¬ 
uid state, and for the alloy during a quench down to 
the amorphous state. Comparison with ab-initio simu- 
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E-Ep[eV] 


FIG. 11. (color online) Instantaneous total DOS in typical 
configurations at 300 K. Dotted line: N = 260; full and dash- 
dotted lines: N = 240. Calculations were made with 36 k- 
points in the tetrahedron method^. 


lation shows that the parameters adjusted in the liquid 
are transferable at high temperature independently of the 
composition, and well below the critical glass transition, 
as far as the pair structure is concerned. This contrasts 
with the parameters adjusted in the solid which cannot 
describe the structure in the liquid, and predict an amor¬ 
phous structure that is inconsistent with ab-initio simu¬ 
lation and the more recent diffraction experiments. This 
rises the question of the path followed to reach possibly 
metastable amorphous states, prior to the comparison 
with theoretical predictions. Concerning the question of 
the size of the simulated system, it is found from a dis¬ 
cussion of the size dependence of the average structure, 
and of the metallicity that ab-initio simulations are fea¬ 
sible on medium size computers while the study of other 
properties such as local structural organization, of low 
cooling rates, or simulations in NPT conditions still re¬ 
quire larger computational resources. Comparison with 
the prediction of the critical temperature of the glass 
transition and the available experimental structure using 
classical simulation, which is less subject to system size 
limitations, finally underlines the importance of develop¬ 
ing thoroughly tested parametrized force fields, in order 
to keep simulations convenient enough for a quantitative 
study of the behavior of complex materials. 
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